COUPLING TECHNIQUES FOR NONLINEAR HYPERBOLIC EQUATIONS. III. 
THE WELL-BALANCED APPROXIMATION OF THICK INTERFACES 
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Abstract. We continue our analysis of the coupling between nonlinear hyperbolic problems across possibly resonant 
interfaces. In the first two parts of this series, we introduced a new framework for coupling problems which is based on the 
so-called thin interface model and uses an augmented formulation and an additional unknown for the interface location; 
this framework has the advantage of avoiding any explicit modeling of the interface structure. In the present paper, we 
pursue our investigation of the augmented formulation and we introduce a new coupling framework which is now based 
on the so-called thick interface model. For scalar nonlinear hyperbolic equations in one space variable, we observe that 
the Cauchy problem is well-posed. Then, our main achievement in the present paper is the design of a new well-balanced 
finite volume scheme which is adapted to the thick interface model, together with a proof of its convergence toward the 
unique entropy solution (for a broad class of nonlinear hyperbolic equations). Due to the presence of a possibly resonant 
interface, the standard technique based on a total variation estimate does not apply, and DiPerna's uniqueness theorem 
must be used. Following a method proposed by Coquel and LeFloch, our proof relies on discrete entropy inequalities for 
the coupling problem and an estimate of the discrete entropy dissipation in the proposed scheme. 

1. Introduction. 

1.1. Main objective. Mathematical models involving a coupling between distinct nonlinear 
hyperbolic systems arise in many applications in physics and engineering science and typically 
involve a non-homogeneous flux function which exhibits discontinuities in the spatial variable. 
Over the past decade, a considerable attention has been paid to the so-called 'conservative cou- 
pling' framework, in which conservation of the unknown is imposed at flux discontinuities 
|HJ |9l [lOj [19J. [20 . 49 1 (and the references therein). Assuming the physical system to be isolated 
in the thermodynamical sense, the conservation requirement is nothing but expected, and many 
problems of interest naturally fall within this framework. In sharp contrast, several other applica- 
tions of equal importance lead to non-isolated systems, which interplay with (possibly singular) 
external sources, the latter (on purpose) locally breaking the conservation property. Typical 
examples are provided by passive or active control devices, while others may fall within this 
category when understood in a broader sense. Considering, for instance, fluid flow problems, we 
observe that momentum and/or energy may be locally supplied or tempered by a wide variety 
of apparatus, ranging from mechanical to electro-magnetical mechanisms. One intends here to 
minimize singular head loss and pressure drop, or accelerate and heat a gas; these apparatus may 
also be used for mixing or cooling purposes. We refer for instance to the book by Gad-El-Hak |31 1 
for a review of current techniques in aerospace and [11] for nuclear safety analysis. Mass may be 
even locally taken from the flow and then injected at a convenient other location [29J in order to 
prevent oil transportation pipelines from slugging. 

In the design of large systems, the fine scale description of the control is commonly bypassed, 
and instead the modeling reproduced its net effect as a sharp transition experienced by the flow 
at the location of the device. The thermodynamical properties of the flow may be affected by 
the control, but even if the flux functions are identical, the resulting jump conditions are not in 
conservation form. Arguments from physics and experiments commonly provide semi-empirical 
laws which express the right-hand trace at the standing transition as a nonlinear function of the 
left-hand trace. This function defines the "transmission conditions" and provides the basis for a 
mathematical formulation in terms of a kinetic function [44 , 45 1 or a family of paths [42 , 26 1 . 
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Various ad hoc numerical methods have been proposed in order to incorporate these trans- 
mission conditions. Recent investigations devoted to a finer design of the transient operating 
conditions of large systems have revealed several shortcomings ITTTI and have assessed the need 
for a rigorous mathematical investigation. As far as the problem of coupling fixed interfaces is 
concerned, the pioneering work [33 j treated the coupling of scalar conservations laws modeling 
the coupling interface as being infinitely thin. Several extensions to the case of hyperbolic systems 
with possibly distinct sizes and involving general transmission conditions have been proposed 
|4 . 15 , 16 1, and transmission conditions are formulated in a weak sense via a 'double' initial value 
problem (IBVP) formalism. They can be as well understood as a measure source term whose 
mass precisely defines the expected departure from conservation, as was proposed in |36| and 
1 44 1 for nonconservative and interface problems, respectively. Various numerical methods have 
been proposed in order to exactly capture (isolated) transmission discontinuities in the setting of 
the coupling problem under consideration in the present paper; these methods are called "well- 
balanced" with respect to the singular transmission source term. (See [2] for related matter and 
11231 for a survey.) 

A difficulty arising with thin interfaces lies in the fact that the initial value problem, even 
with apparently well-defined interface conditions, may turn out to be ill-posed, so that the thin 
interface model does not fully determine the dynamics of the fluid flow. This feature is related 
to the resonance that may take place at the interface, when waves associated with the fluid have 
almost vanishing speed. Even in the scalar case 1141 , multiple solutions to the initial value problem 
are exhibited when the coupling interface is resonant. The failure of uniqueness corresponds in 
fact to a general situation first described in [37 1 and further analyzed in [32J. 

The present work is a continuation of our analysis in 1 1611171 (to be continued in [18 1), which is 
devoted to resonant coupling interfaces. In the first two parts of this series, we introduced a new 
framework for the mathematical coupling, based on an augmented formulation which has the 
advantage of avoiding any explicit description of the interface structure and was referred to as the 
thin interface model. The coupling problem takes the form of a standard IBVP problem, which 
can in turn support various regularizing mechanisms. In [16 1, we relied the self-similar viscosity 
method by Daf ermos 1 25 1 and established the existence of self -similar solutions (with shock waves) 
for the coupling problem of two hyperbolic systems (under fairly general assumptions). However, 
in the limit of vanishing viscosity parameters that we studied in Ifl7l , a lack of uniqueness 
is observed for solutions involving a resonance effect, even in the simple scalar setting. We 
emphasize that entropy inequalities that would attempt to incorporate at the macroscopic level 
the fine scale effects modeled by viscous mechanisms, do not lead to a efficient selection principle 
for thin coupling problems. 

In the present paper, the augmented approach proposed by the authors [ 16 1 is shown to lead 
to another regularization strategy, now based on a thick interface model, as we call it. Roughly 
speaking, the singular source term modeling the transition is given a smooth profile but the 
overall regularization technique achieves the key property that the left- and right-hand traces of 
any isolated transition waves are still exactly captured. Since the source term is entirely localized 
within the transition wave, it does not act elsewhere and steady solutions of the IBVP problem 
are thus expected to stay constant outside the transition profile. This assesses the importance 
of considering isolated transition waves. Importantly, this accuracy property holds for resonant 
transition waves. It is achieved thanks to a well-balanced strategy that can be traced back to the 
seminal work by Greenberg and Leroux [35 1 (see also Bouchut [13 1, Gosse [34] and the references 
therein). This well-balanced property holds for any given regularized profile, in the setting of 
general pair of fluxes and transmission conditions. The versatility of the method allows us to 
address a fairly general non-conservative coupling problem for which a given smooth profile 
may be promoted from experiments and knowledge considerations. 

An outline of this paper is as follows. Focusing on equations in one space variable, we 
briefly recall some of the existing frameworks for the non-conservative coupling of two scalar 
laws. We then introduce the augmented PDE model with thick coupling interface. Existence and 
uniqueness of an entropy weak solution follows from the well-known Kruzkov's theorem. We 
propose a scheme for approximating the solution of the Cauchy problem and prove the expected 
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well-balanced property. Due to a reconstruction procedure of the discrete solutions dictated by 
the well-balanced property a priori uniform BV (bounded variation) estimates are not known in 
general and strong convergence is proved by following the entropy dissipation method of Coquel 
and LeFloch [24], using DiPerna's uniqueness theorem in the class of entropy measure-valued 
solutions. Finally several numerical results illustrate the flexibility of our strategy and its ability 
to capture various resonant transition waves. 

1.2. Background. We consider the coupling between two conservation laws at an interface 
located at the location x = on the real line R: 

(1.1) d t w + d x f ± {w)=0, ±x>0, 

with unknown w = w(t,x) e ]R, defined for t > and x e ]R, where the fluxes /* are twice 
differentiable functions. Prescribing the initial data iv(0,x) = Wq(x) at time t = does not suffice, 
and an additional condition for modeling the transient exchange of informations at x = must be 
supplemented. Such an additional closure is called a coupling condition and is motivated by the 
general considerations made in the introduction. Focusing at this stage on piecewise Lipschitz- 
continuous solutions with bounded left- and right-hand traces, this coupling condition can be 
given the form of a nonlinear closure law expressing on the left-hand state as a function of the 
right-hand state of vice-versa. Many different coupling conditions may be introduced and should 
reflect the precise departure from the conservation property observed in specific applications. 
After |3|-|7J and [22J, the coupling condition is formulated from two monotone (increasing, say) 
functions ± :ib€Kh 6±{w) 6 M, so that one writes 

(1.2) (0_oH?)(f,O-) = {6 + ow){t,0+), f>0. 

Two approaches are available for general coupling problems, as now discussed. 

On one hand, infinitely thin interface models are based on prescribing coupled boundary 
conditions, as pioneered by Godlewski and Raviart [33J; see also Boutin et al. [14J. Boundary 
conditions are formulated so that dL2j is understood in a weak form, following Dubois and 
LeFloch 1 28]. In order to formulate the Cauchy problem in the half-space M + , a boundary 
condition feiteRn b(t) e ]R is prescribed at x = and imposed in the sense w(t,0+) e + (b(t)) r 
where 0+(b(t)) = |W(0+, b(t),w),w e mJ denotes the set of admissible traces W(0+, b(t),w) at 
x = 0+ of Riemann solutions W(-,b(t),w) associated with prescribed left-hand state b(t) and 
arbitrary right-hand state. Denoting by 0+ the set of admissible traces at x = 0± determined from 
Riemann s olut ions with flux f* and boundary states o0 ± o w(t, 0±), we formulate the coupling 
condition |l.2| in the weak sense w(t,0±) e Q±(dl 1 ° 6± o w(t,0±)j,t > 0. It can be checked flit ) 
that the proposed weak form for the coupling condition reduces to the strong version | |1.2| as 
long as a resonance phenomena does not take place at the interface. By resonance, it is meant that 
waves from either the left- and/or right-hand problems interact with the interface, so that the 
continuity property | |1.2| is lost in general and multiple discontinuous solutions may be available. 
Therefore, an additional selection criterion is required. 

1.3. Augmented model with thin interface. On the other hand, in the second framework 
for coupling problems, the thick interface model of the authors 1 16 1 (also considered earlier in 
1 15 1 for the self-similar regularization of scalar equations), we view the coupling interface as a 
standing wave for an augmented system of partial differential equations, by generalizing here 
the nonconservative reformulation proposed by LeFloch IBTI [42] for the nozzle flow problem 
with discontinuous cross section (cf. also [46 J). The standing wave is designed so that (away 
from resonance, at least) a complete set of Riemann invariants is available, in agreement with the 
continuity property ^1-2} . We propose to consider a new unknown u = u(t,x) defined by (t > 0) 



(1.3) u{t,x) 



(d-ow)(t,x), X < 0, 
{6+ow){t,x), x > 0, 



so that the strong coupling condition 1 1.2 1 is equivalent to the continuity condition 
(1.4) u(f,0-) = M(f,0+), t>0. 
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Let us stress from now on that this convenient reformulation of the coupling condition will play 
a central role in the derivation of a well-balanced method. Observe that < |1.3) is a well-defined 
change of variable, since d w 6 ± {w) > for all aieR, We introduce an unknown v = v(t,x) which 
coincides with the Heaviside function for all f > 0: 

(1.5) v(t,x) = vo(x) = if x < 0, 1 if x > 0. 

The value is meant to recover the given equation in the half line JR~, while the value 1 represents 
the given equation in ]R + . Intermediate values of V model a smooth transition region from one 
problem to the other. We thus introduce the augmented model (f > 0, X e K) and its associated 
initial condition 

d t C (u, v) + d x Ci(u, v) - d v Gi(u, v)d x V = 0, d t V - 0, 

(1.6) M _ J (0- o HToKac), x < 0, _ / 0, x < 0, 
" oW " \ (0 + o Wo )(x), x > 0, y ° (x) " \ 1, x > 0. 

Here, denotes the initial data for the coupled problem | |1.1) while the functions Go and G\ are 

(1.7) e (u,v) = (l-v)y-(u) + vy + (u), Qi{u,v) = (1 - v)f-(y.(u)) + vf + (y + (u)), 

and y± are defined to be the inverse functions of the increasing map 8±, respectively. By our 
monotonicity assumption on 8±, one has 

(1.8) d u Q {u,v) > 0, ueK, ve[0,i\- 

This property obviously preserves the time direction determined by the nonlinear first-order 
augmented system | |1.6) . Other choices of the coupling functions So, Ci are possible, as we 
discussed in lTl6l . Observe that smooth solutions to l |1.6| obey 

(1.9) d u e (u,v)d t u + d u Qi(u,v)d x u = 0, d t v = 0, 

so that the first-order system admits two real eigenvalues: and A(w, v) = (d u Go(u, v)j d u Q\{u, v). 
This system also admits a basis of eigenvectors and the characteristic field associated with A(u, v) 
is genuinely nonlinear, provided the flux functions / ± are genuinely nonlinear. The other field is 
linearly degenerate and the standing wave is clearly characterized by d u Q\(u,v) d x u = 0, which, 
provided that d u Ci(u, v) + 0, implies the Riemann invariant u = est, so that the coupling condition 
m(0-, t) = u(0+, t) stated in |L2| is satisfied in a strong sense. 

States (u+fV*) with d u Q\{u* ,v*) = may exist when one (or both) speed (/*)' changes sign. 
Observe that such states come with the property A(u*, v*) — and thus correspond to the interac- 
tion of a possibly genuinely nonlinear field with a linearly degenerate one. We must then define 
weak solutions to the non-conservative nonlinear system | |1.9) when eigenvalues vanish. Solutions 
for such hyperbolic systems are not uniquely defined, unless additional physics is prescribed, as 
recognized in LeFloch lHTll42l|45l . This situation is referred hereafter as to a resonance phenom- 
ena. Indeed it has direct connection with the setting for resonance investigated by Isaacson and 
Temple [37J and Goatin and LeFloch Il32l . As already reported, the definition of weak solutions 
for | |1.6[ in the resonant regime has been tackled in [15] via the Dafermos self-similar vanishing 
viscosity analysis. In |17J, multiplicity of self-similar solutions is shown to persist for the non- 
conservative model in the limit of vanishing viscosity. Failure of uniqueness arises for resonance 
problem as noted in |37| and for various interface problems even in linear hyperbolic equations 
|45 Chap. 5]. The origin for multiple solutions is found in the property that Riemann solutions 
describe the time-asymptotic behavior of the Cauchy problem for parabolic perturbations of | |1.9[ . 
Here, the precise definition of a regularization v n (for some ry > 0) plays a central role in the non- 
uniqueness of Riemann solutions. The regularized profile v J] does not weight the wave speeds / ±/ 

equally within the expression \{u,v n ) = (f? i( eo(«,c, ; )) ((l-v 1l )y'_(u)f~'(y-(u))+v, l y' + (u)f + '(y + (u)fj. 
Consequently, in the resonance phenomena more importance is given to the left or right-hand 
problem and this is the origin of the failure of uniqueness. We refer the reader to |Tl7l where up 
to four solutions can be build from self-similar analysis. We provide below numerical evidences 
that multiple solutions are stable, in the sense that each solution can be captured numerically. 
(See also Schecter et al. |47. 48] for a discussion of multiple self-similar solutions.) 



1.4. Thick interface model. We now extend the previous step-like color-function to a smooth 
color-function v = v{x) and we consider 



(1.10) 



w(t,x) = e (u(t,x),v(xfj, t>0, xeK. 



In view of the monotonicity property | |1.8) satisfied by Co(., i>)> the function u can be recovered 
from w, the color function v being fixed. With some abuse in the notation, we write w — w(u, v) 
and u = u(w, v). The interest in this change of variable stems from the fact that the first equation 
in \l-9\ reduces to the balance law describing the thick interface model 

w(0,x) = IVq(x), 



(1.11) d t w + d x f(w,v) = £(w,v)d x v, 

where, in agreement with 1 |1.6) , 
(1.12) 



f(w,v) = Gi(u(w,v),v), £{w,v) = d v Q\{u{w,v),v) 

and the initial data is w = w(uq,v), in agreement with | |1,6} . Let us stress that the product 
t(w, v)d x v is now nothing but a standard zero-order source term, due to the smoothness of v. 

The thick interface framework for coupling problems allows us to use the notion of entropy 
pairs for conservation laws with source terms. Any convex function XL = U(w) can be used 
to define an entropy pair (II, GF). To define the required flux, we start from the augmented 
system | |1.6) and write, for smooth solutions, dt£o(u, v(x)) + d u Gi(u, v(x))d x U - 0, which leads us to 
d t U(e Q (u,v)) + d x Q(u,v) - d v Q(u,v)d x v = 0, with 

(1.13) Q{u,v) = J U'(e (6,v))d u e 1 (d,v)d6. 

In terms of the unknown w, this reads d t U(w) + d x 3^(w, v) - L{w, v)d x v = 0, with 



(1.14) 



3(w,v) = Q(u(w,v),v), 



£>{w,v) = d v Q(u(w,v),v). 



Weak solutions of the conservation law with smooth spatial rnhomogeneities 1 1.11} are then 
naturally selected by the entropy inequalities 



(1.15) 



d t U{w) + d x 3{w,v) - L(w,v)d x v < 0, 



understood in the distributional sense for any convex entropy pair (IX, J). Here and since again v 
is smooth, L(w,v)d x v acts as a usual source term and Kruzkov's theory [40 1 applies and provides 
us with a unique entropy solution to 1 1.15) , when the flux and source terms are sufficiently 
regular, say piecewise differentiable. The minimal smoothness property on the color function 
v to meet the Kruzkov assumptions is therefore v e W 2oo (]R). We will see that existence and 
uniqueness of an entropy solution of I l.ll) - jT7l5| in fact holds (under this smoothness condition 
but) for general initial data, that is, Wq e L°°(R). 

In order to motivate our method, recall here some properties of time-independent solutions to 
1 1.11 1, i.e. solutions satisfying d x f(w, v) = t(w, v) d x v or, in the w-variable, ((1 - v)f~' (y-(u))y'_(u) + 
vf + (y+(u))y' + (u))d x u - 0. At the numerical level, it is very challenging to capture steady solutions, 
especially when the coefficient ((1 - v)f~'(y-)y'_(u) + vf + '(y + )y' + (u)) vanishes, that is, when the 
non-trivial eigenvalue A(u,v) of the hyperbolic system 1 1.6 1 vanishes — which is the resonance 
phenomena. In view of the coupling condition 1 1.4 1, our strategy will be to focus on (non constant) 
solutions w to 1 1.11 1-1 1.15 1 which have constant component u but variable component v which 
are treated as "stable solutions" — even when resonance occurs. 

In the next section, we therefore introduce a well-balanced finite volume method for approx- 
imating the entropy solution to | |1.11 l- |i~15 1. As already stressed, by well-balanced we mean 
that solutions w{u,v) with constant components u and general components v, so that u{w, v) is 
constant in space and time. An adapted reconstruction procedure will be required in order to 
achieve our goal. The resulting family of approximate solutions will be seen to be uniformly 
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bounded in sup-norm under a natural CFL (Courant-Friedrichs-Levy) condition. However a 
uniform a priori estimate in total variation is not available (except for the trivial case w(u, v) = u), 
due to the well-balanced reconstruction step. Following Coquel and LeFloch E4l , we propose to 
take advantage of the existence of infinitely many entropy differential inequalities and advocate 
the use of DiPerna's theory of entropy measure-valued solutions ||27| . 

2. Formulation of the well-balanced scheme and convergence theorem. 

2.1. Notation and assumptions. In this section, we present a finite volume method for the 
approximation of entropy solutions satisfying, by definition, and ^1.15} . As explained, 



we require the method to be well-balanced with respect to the family of stationary solutions 
w = zv(x) relevant for the coupling problem, that is, solutions characterized by the condition that 
u{w{x),v{x)) is constant in x e R. Since v depends on x, so does the stationary solutions w{x). To 
ensure the well-balanced property, it is convenient to design a finite volume method and handle 
the two components of the system on two distinct grids. (For this strategy, we refer the reader to 
the review (T3|, as well as l2Tll30ll39ll). 

For simplicity and without genuine loss of generality, we consider constant time and space 
steps, denoted by At, Ax > 0, respectively. We then introduce the time levels t" = nAt (n = 0, 1, . . .), 
the cell centers Xj = jAx, and the cell interfaces Xj+i/2 = (j + 1/2) Ax (for all integers /). The 
approximate solutions u Ax and v Ax are sought as piecewise constant functions, with 

u Ax (t, x) = m* x e {Xj- VZ , Xj-yz), t e (t n , t n+1 ), 

VAx(t,X) = Vj+i/2, xe(xj,x j+1 ), f>0, 



and, in view of ( 1.10} , we have also the companion function 



(2.2) w Ax (t,x) = w(u Ax (t,x),v Ax (t,x)), t > 0, x e R. 



Since the solution v in ^1.6} is independent of the time variable, v Ax (t,x) is chosen to be time- 
independent. We also set 



I rXj+i/2 

M L( X ) = u °j ~ t - I u Q (x)dx, x e (Xj-y 2 ,Xj- 1/2 ), 

V Ax {t,X) = Vj+l/2 = — J V(x)dx, X e (Xj,Xj +1 ). 



(2.3) —-v.,- i; 



The discrete solution u Ax will be evolved in time by a finite volume method, consistent with the 
equation of interest 

(2.4) d t w{u,v) + d x f(w(u,v),v) - (f + {y + {u))-f-{y-{u)))d x v = 0. 

By construction, the discrete function v Ax is constant within a neighborhood of each cell interface 
*7+l/2/ so ma t the above equation locally reduces to a conservation law in the unknown w = 

w(u,Vj +1/2 )- 

(2.5) d,w + d x f{w,v j+vi ) = 0, xe{xj,Xj+i), te(t",t" +1 ). 

This property motivates us to introduce, at each cell interface Xi+yz, a two-point numerical flux 
g(-, •; C/+1/2) : R x R -> R, which is assumed to be locally Lipschitz continuous and satisfy the 
consistency and monotonicity properties 

(2.6) g{a,a;Vj +V2 ) = f{a,Vj+y 2 ), a e R 



(2.7) j-{a,b;vj+y z ) > 0, -±{a,b;v j+ y 2 ) < 0, a,b € R. 
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2.2. The well-balanced scheme. The discrete solution u&x(t n , .) being known at time f", we 
determine the new approximation at the time t n+1 into two steps: a subcell reconstruction step 
followed by an evoJution step. In turn, our algorithm is a time-explicit finite volume method, 
which will be shown to converge under the (CFL) stability condition 



(2.8) 



At 

— max sup 

A* i ue[m,M] 



df 

— (w(u,V j+1/2 ),V j+V2 ) 



1 

*2' 



where m = inf YSEE u (x) and M = sup x€K u (x). 

• Subcell reconstruction. At each time t" and in each cell (xj-i/2,Xj+\/2), using the change of 
variable w - w(u, v) we introduce two "subcell states" v^Ly 2 ± to g e t ner with their average: 

(2.9) w n hllZi+ = w(u'j, v hV2 ), w n j+vx _ = w{u n r v j+1/2 ), w'j = l(w1_ 1/2+ + w n j+V2i _\ 

• Evolution in time. At the time t n+l and in each cell (xj-i/ 2r Xj + i/ 2 ), we define w". +1 by 
integration on subcells and set 



(2.10) 



io n+1 = w" 



(G» 



7+1/2,- 



With G " + l/2, 



S( w1 j+i/2,-' wr j+i/2,+' v i +l/2 ^ ~ f( m 1+V2-' v 1+V 2 ) and a similar expression for 
G"-i/ 2 +■ Then, the "new state" u" +1 is defined as the solution to the algebraic equation 



(2.11) 



\(w{uf\v hll2 ) + w{uf\v M i 2 )) = wf\ 



This completes the description of the proposed method. 



The monotonicity property | |1.8) , namely d u w(u,v) > 0, ensures that (2.11 1 admits a unique 
solution. Observe that < 2.10) is also equivalent to 



with gy +1 / 2 ~ g( w 'j+i/ 2 -' w "j+i/2 +' p /+l/2) we obtain from the definition 1 1.12 1 of f(w, v) 

f{w" +ll2i _,v j+V2 ) - fW\- V 2,+> v i-vi) = (/+(y+(w")) - f-{y-{u"))){vj + in - v hV2 ), 

so that 1 2.10| is a formally consistent discretization of the governing equation 1 2.4 1 in the unknown 
w. 

2.3. Main convergence and well-balanced results. We can now state our main results. 

Theorem 2.1 (Convergence of the finite volume method). Consider the Cauchy problem ( |L6| 
with initial data uq in L°°(]R) and vq in W 2 '°°(R, [0, 1]) and suppose that the monotonicity property 
d u w(u,v) > holds (see (1.8 Consider the family of functions v Ax defined in (2.1 H2.3 1 and the family 
of approximate solutions w^ x defined by \2.2\ (from u^ x in \2.1\ ) and the finite volume method (2.9 H |2 JO 1 
whose numerical functions satisfy §2.6\ -(2.7K Then, under the CFL condition (2.8* and as Ax —> 0, the 
solutions if a x remain bounded in L°°(]R + x R) and converge strongly in the LfLT norm (1 < p < oo) 
toward the unique Kruzkov solution w to (1.11 *- \1.15 1. 

Proposition 2.2 (Well-balanced property). Consider the Cauchy problem (1.6 1 when the the initial 
data Uq(x) = u+ (x e is a constant and the data vo : R — > [0,1] is any smooth function. Then, the 
discrete solution u^ x given by \2.9\ -(2.10* is also constant in space, with u& x (t n , x) = u+ (x e Mj at each 
time level t n . 

Proof. The discrete initial data is now u°. = u+ for all j so, at each interface Xj+i/ 2 , we 



get w 



7+1/2,- 



- W° 

;+i/2,+ 



= w(u+, Vj+1/2), irrespective of the state Vj+\/2- The subscript ± may be 
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omitted and the flux at Xj+1/2 reads g° 



/+1/2 



(w° j+1/2f w° j+V2 



;v 7+1/2) 



7+1/2 



/ £7+1/2) (in view of 



the consistency property |2.6l). As a consequence, the left- and right-hand fluxes 12.10 1 vanish 
identically: GJ +1 ^ 2 _ = G" +1/2+ = for all /'. The scheme 12.10 1 thus yields w l . = zv°. for all ; and, in 

view of 1 2.9 >, u 1 . in each cell satisfies io{u l -,Vj-\iz) + w(uJ/uJ^i/ 2 ) = w(u+,Vj-\i 2 ) + if(i/*,c /+ i/2). In 
view of the monotonicity property 1 1.8 1, namely d u w > 0, we deduce u l . - u+, and an induction 
yields the desired conclusion. □ 

2.4. Formulation based on convex combinations. We briefly revisit the finite volume method 
( |2.9) -| 2TT0) so as to highlight its relationships with existing well-balanced approaches. Then, we 
put forward a convex combination (at the level of subcell) which is of central importance in our 
forthcoming analysis. Consider the following auxilliary Cauchy problem (t e (0, At), x e R) 



(2.12) 



d t w(u, v) + d x f(w(u, v), v) - €{w{u, v), v)d x v = 0, 



d t v = 0, 



(2.13) 



(u Q (x),v (x)j = [u hx {f,x),v hx {x)) 



_j (u",Vj- 1/2 ), X £ (x h i/2,Xj), 

(u",v j+1/2 ), x e (xj,Xj +V2 ). 



Solving this Cauchy problem in the time slab (0, At), with At satisfying the CFL restriction l |2.8} , 
just amounts to glue together non-interacting Riemann solutions emanating from the interfaces 
Xj and xy+i/2. Indeed, at Xj, the Riemann data has an arbitrary jump in Vq(x) at x - but with 
a constant u${x) = u'j, this property allows us to solve the (local) Riemann problem in term of 
a standing wave (depicted in Figure 2.1 as a vertical line emanating from Xj). Here, we thus 
favor such a stationary solution even if resonance locally takes place. On the other hand, at xy+i/2 
where vq(x) is locally constant, the Riemann solution u is easily determined at time At, that is, 
w(u(At,x),Vj +1/2 ) = M (tt^ w " +1 / 2 - r w " +1 / 2r+ ) where <y(f; w" +1/2 w" +1/2 + ) is the self-similar entropy 
solution of 



(2.14) 



d t co + d x f(a>,Vj+i/2) = 0, t > 0, x e R, 

' W 'j+l/ 2 - = w ( u ",Vj+l/2), X < 0, 
W ;+l/2,+ = W ( U "+l> V i+V2), X>0, 



co(0, x) = zv (x) 



in agreement with the subcell reconstruction step 1 2.9 1. Observe that under the CFL condition 
|j2.8}, the solution u(t, x) of 1 2.12^ cannot interact with the two neighboring standing waves at Xj 
and Xj+\ for times t e (0,At). Hence the exact solution to the Cauchy problem 1 2.12 1-12.13 1 is 
obtained in the desired form (Cf. Fig. 2.1 1. 




projection 



advection 



reconstruction 



Figure 2.1: A subcell convex combination 
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Next, note that an averaging procedure at time At in each subcell (Xj-\/2,xA and (xy, x ;+ i/2) 
yields 



w 



,.n+l- 
/-l/2,+ 



2 

Ax 



r x i 

I w(u(At,x),Vj- 

•JXj-lll 



1/2 



)dx, 



,.n+l,- 
7+1/2,- 



2 

Ax 



1 1 



w(u(At, x), Vj+i/2)dx. 



Using classical arguments, the arithmetic average of the above subcell states gives 

At : 



(2.15) 



where in agreement with (2.10 -1 2.11} , the left- and right-hand Godunov fluxes read GJ +1 ^ 2 _ = 

/(fl)(0+; Wj +1/2 IP ■ +1/2/+ ), ^W-/T^i/2,-' ^7+1/2)' 311(1 G *-l/2,+ = ^-1/2,-' W "-l/2, + )' ^7-1/2)- 

f( w 'j-i/2+' v j-~t/2)- m other words, the formula (12.10 for w" +l is recovered in the special case of 
the Godunov solver for system (12. 12b . Other lett-hand and rieht-hand fluxes (based on gen- 



eral monotone numerical flux functions g(., .;cy+i/2) satisfying fl2.6} -|2.7 <) are also obtained by 
approximating the Riemann solution ( 2.14} . We summarize our result as follows. 

Lemma 2.3 (A subcell convex combination). Consider the finite volume method (2.9 H2.20 1 with 
fluxes g(., .;Vi+uz) satisfying ( 2.6. H 2.7,1, and introduce the subcell states 



(2.16) 



w 



,.n+l- 
7+1/2,- 
n+1,- 



W ;+l/2,- 



^Vf/2,+ = W /+l/2,+ 



2Af 
Ax 
2Af 
Ax 



(^+1/2-/(^+1/2,-^7+1/2))' 

r(/K+i/2,+'^i/2)-^ + i/ 2 )' 



ipz/-/i g\y 2 = g{ w " + \i2 -> w 'j+\/2 +' /V 7+i/ 2 )- the formula {2.10 \ for w recasts in term of the convex 

combination 



(2.17) 



1 2 V "7-1/2,+ /+1/2 -J- 



Rephrasing the above result, the formula (2.10 1 for zv" +1 , which we have underlined to yield 
a consistent discretization of the equation 1 2.4} , just reads as a convex combination but for the 
conservation law | |2.5} . The next equality reflects such a property and it will be extensively used 
in the sequel 



(2.18) 



11+1 



= \{w{uf\v hl i 2 ) + w{uf\v j+l i 2 )) = l(wp /: 



1/2,+ 



■ W 



n+1,- 
7+1/2, 



3. Well-posedness theory for the thick interface model. To motivate the forthcoming de- 
velopment, we first recall that in the present coupling setting with thick interfaces, the color 
function v is a given smooth function of the space variable, say v e W 2,00 (R). Hence as already 
emphasized, the Kruzkov's theory applies to establish uniqueness of the entropy solution w in 



L™ (R+ x R) of the Cauchy problem |TTTJ-([Tl5} with initial data iv e L°°(R). We will establish 



hereafter the required sup-norm estimate for the family of approximate solutions Wax, with Wax 
defined in (12.2b, but no such uniform estimate is known in the BV semi-norm. Indeed, the total 



variation of the discrete solutions may increase at the subcell reconstruction step (2.9 1 (except in 
the particular case y + = y~ = Id, see section |5| and a control of the total variation seems out of 
reach. The absence of an a priori strong compactness argument leads us to adopt the setting of 
measured-valued solutions for ( l.ll} -l i~T5} so as to recover, following DiPerna ||27|| , a posteriori 
strong convergence of the approximate solutions w& x on the basis of infinitely many entropy 
inequalities. 

In this section, we state a generalization of DiPerna's uniqueness theorem 1271 which concerns 
the class of entropy measure-valued solutions to nonlinear equations | |1.9) . Measure-valued 
solutions are Young measures, that is, weakly measurable maps /./ : (f, x) e R+ xR-> ji t ,x which 
take their values in the space of probability measures and, in our case, are supported in a compact 



interval of R. A Young measure represent all weak-star limits a(w^ x ) of a bounded sequence w& x 
for arbitrary functions a € C°(R), that is, a(w^ x ) — 1 (/.; t v ,a) = a(A)d/i (x (A) weakly - ★ in L°°. 

Definition 3.1 (Entropy measure-valued solutions). Let v 6 W 2/ °°(R) and u e L°°(R) be given, 
and let wq = w(uq, v) be an initial data for the Cauchy problem {1.11 K A measure-valued map \i tiX is called 
an entropy measure-valued solution to the Cauchy problem \1.11\ if for every convex entropy pair (U, IF) 
in the form {1.13^-(l.li\ one has 

f ((Htji, U)d t (p + {y. tlX , v))d x <p\ dtdx + I cp(^, £»(; v))d x v dtdx + I U(w )(p(0, .)dx > 

Jr+xR V ' Jr+xR Jr 



Rwe set {nt A ,X{',v)) = 



for any non-negative test-function <p e f(]R+ x R). 

With obvious notation, for any continuous function J£ : 1R X [0, 1] 

Theorem 3.2 (Uniqueness in the class of entropy measure-valued solutions). Let v be giv en in 
W 2 '°°(R), uq in L°°(R) and fi = y. t/X be an entropy measure-valued solution (in the sense of definition 3.1 1 
of the Cauchy problem \1.11 1 with initial data wq = w(uo, v). Then, for almost every (t, x), the measure [it, x 
is a Dime mass pL tiX = 5 w (t, x ), where the function w e LT (R+ x R) denotes the unique Kruzkov's solution 
of the Cauchy problem \1.11 H 1.15>. 

The proof of this result follows from the one in Ben-Artzi and LeFloch in [12 j and Amorim, 
LeFoch, and Okutmustur |8J for conservation laws on manifolds and details are left to the reader. 
Observe here that the initial data is automatically assumed in a strong sense, namely 



(3.1) 



lim I I (u t x , |Id - w \) dtdx = 
!?° Jo Jk 



for all compact K in R. The following technical lemma provides us with this property, while 
f further material can be found in ft2"ll50ll . 

Lemma 3.3 (DiPerna [27]). Suppose that there exists a strictly convex function U and a Young 
measure f.i satisying, for allip >0 in C£°(]R), 



(3.2) lim I I ([i tiX ,ld)il> dtdx = I wo^>dx, lim I I (j.t t , x ,U)ip(x)dtdx < I 
"o Jo Jr Jr Jo Jr Jr 



U(wo)ipdx. 



Then, the property {3.1 1 holds. 
4. Convergence analysis. 

4.1. Local maximum principle. We show now that the approximate solutions remain bounded 
inL°°(R + xR). 



Proposition 4.1 (Uniform sup-norm stability). Under the CFL condition \2.8\ , the finite volume 
method {2.9 H2.10 1 satisfies 



(4.1) 



min( M J_ 1/M ^Mj +1 ) < uf' < max{u n hV u n .,u n j+l ) 



at all time level t" and, consequently, ||waxIIl»(r + xr) < II"oIIl~ (r) an d ||waJIl~(r + xr) < 0(1). 

Proof. From the subcell states introduced in (2.16 \, let us define the auxiliary quantity 



^"+1/2 - as me unique solution of w 

.n+l 



n+l,- 
7+1/2- 



= w(u"^y 2 _,Vj+i/ 2 ) as well as u 



n+l,- 



w ( u j+l'/2 +' V j+V2)- With these definitions, the identity 12.18 gives w(u" 

n+l 



= w(u n +l'~ + , Vj-nz) - w(u'l +1 ,v hl/2 )- We thus deduce (u? +1 



monotonicity property satisfied by w(., v), that is 
(4.2) 



u n+l 
M ;+i/2 



n+l,- 

/+l/2,+ """" - ;+l/2,+ 
" +1 ^l/2)-W(u^_,0; + l/2) 



from zt> 

/-| 

•7+1/2,-' °/+i/2) 
mJ +1 ) > from the 



min( M ^- + ,„^-_)< M »- 



<max( U ^- + , U »^-_). 
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It suffices to check that the states u "+lj 2 ± satisfy min(«", u" +1 ) < u"*y 2 ± < min(u" ,u" +1 ). For mono 



«+!,- 



tone schemes and under the CFL restriction ( 2.8 >, both subcell states a;"^ 1 ' in (2.16 1 satisfy 



/+l/2,± 



(4.3) 



min(^ 1/2 ,_,^ +1/2j+ ) < w"ty 2A < ^(P"+i/2,-' w hnA 



Thanks to the monotonicity of u(.,v) 
< u(w n+1 - 

and u n - +l = u(w^ +1/z+r Vj + i /2 ) in view of |Z9l. 



w 1 (.,v) (for o 6 [0,1]), we deduce that min(«J, 



; j 2 ± , C/+1/2) < max(«J, «" +1 ), since, in the reconstruction step, one has u" = u(w" +1/2 



□ 



1 are deter- 



4.2. Discrete version of the entropy inequalities. Since the subcell states w'. +l , 
mined from monotone fluxes, they satisfy a discrete version of the entropy inequalities coming 
with ( 2.5 1 and we rewrite here d t w + d x f{w,Vj+\i2) = 0. This is the matter of the next statement 



which is of central importance. 

Lemma 4.2. Let U(.), 3(.,Vj+y2) '■ R —* RxKte any convex entropy pair for the conservation 



law (2.5). Then there exists a locally Lipschitz continuous, entropy flux S(.,.;i> 1+1/2) : RxR -> R 
satisfying the consistency property 3(a,a;v '.+1/2) = 9(a,Vj+i/2) ( a e R)/ so that the following discrete 

n+l- 
j+l/Z,± 



entropy inequalities hold (under the condition (2.8 > and with w"^' n ^ introduced in §2.16}) 



(4.4) 



U(w»;lf 2 j - U(w» +V2 _) +2£(s? +1/2 - J(zv'; +V2 _,v l+V 2)) < 0, 
~ U(w1 +1/2+ ) + 2$±(j(w1 +y2+ ,v j+ y2) - SJ +1/2 ) < 0, 



This result is classical and the proof is omitted. The proposed inequalities are the cornerstones 
of the following (discrete in time but continuous in space) entropy inequalities. 



Lemma 4.3 (Shifted time discrete entropy inequalities). For every test function cp > in CD(R+ x 

R), one defines (p" +1 / 2 = AtKx ft» Ix' + ^ <p(t> x ) dtdx. Then, under the condition (2.8 1, the (discrete-in-time 
and continuous-in-space) inequality 



i i 
( 4 - 5 ) 'J J (^ W lx> v ( x )) d x4>(t, x) + (p{t, x)L(w n Ax , v(x))d x vj dtdx 

< 0(Ax)At\\cp\\ w t^ m n*i )xn) x<p(t n , 
holds, where x$(t n , t n+1 ) = 1 if max t€ ^„ t ,:+^(max x€ K (p(t,x)) + 0, while X<t>{t n , t n+l ) = otherwise. 



4.1 



is 



In addition to the property v e W 2 '°°(R), only the sup-norm estimate in Proposition 
needed to deduce | |4.5) . The latter can be handled in the limit Ax — > 0, by using the Young measure 
H associated with w^ x . To evaluate the discrete time derivative in | |4.5| , we introduce a cell average 
representation of (p, i.e. 



(4.6) 



<P] = l(<Phn + <P} + V2)> 
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and we recall 

£ \(^Z t ,-W h XI2 + "K- + l/2, + )^/- 1/2 )^ " £ \(^ n j+V 2,W i+ll 2 + U(^_ 1/2 , + )^_ 1/2 )A* 

= £ ((u(^ ,+i ) - u(07p)^Aac - \Y, (2 - t^ + i; 2 j^ +1/2 - U(^- + )^_ 1/2 )A^ 

/ i 

+ 1 £ (2 U(^ - U(^ +1/2i _)^ 1/2 - U(^_ 1/2i+ )^ 1/2 )ax. 

(4.7) 

The first term in the right-hand side of 1 4. 7 1 yields the time derivative, and we control the 
remaining term as follows. 

Lemma 4.4. The following estimate holds in each mesh cell 

(4.8) 2 U(w^ - U(w1 +V2i _)^ +V2 - U(w1_ mi+ )^_ V2 < 0((Ax) 2 )\\<p\\ m 

while 

\\dx(p\\L-»((t",t<^)x(x H/2 ,x j+1/2 )), 

where o u denotes a convexity modulus ofU : U"(w) > o u > for all w such that \w\ < Ww^Wico^+^y 

With a strictly convex entropy having ou > 0, the bound < |4.9) is slightly sharper than the 
estimate 



2 U(w^W - U{w^-_W ]+V2 - U{w n ^- + )^_ V2 



(4.10) 2 U(it/' +1 )<^ - U(w^'- 2 Jcp» +in - U(w n ^l + )cpi_ V2 < 0(Ax)\w n ^f z _ - w n ^ + \ ||^|| L ~, 

which is a crucial observation. The derivation of ( |4.9) relies on the following result whose proof 
is left to the reader. 

Lemma 4.5. Let w : x e (0, 1) — » zi>(x) e ]R be a bounded function with mean value w. Then, for any 
convex entropy pair U:!»eR-» U(a;) e ]R the following estimates hold with m = min A - e (o,i) w(x) and 
M = max xe(0 ,i) w(x): 



(4.11) min 

m<v<M 



w(v) r 1 r 1 

— ^ \zv(x)-W\ 2 dx< 
2 Jo Jo 



lX(a;(x))dx - ll(it>) 



U"(i>) f 1 
< max — - — I |zt>(x) - w\ dx. 

m<v<M 2 Jq 



Let us comment about | |4.8) and | |4.9) in Lemma 4.4 with a somehow sharper version of ( |4.8) 
one has 

2 U(w^» - U{w" j+V2 _)^ +ll2 - U(w1_ V2)+ )<p»_ y2 < 0(Ax)\w'] +1/z _ - w n hV2i+ \ ||^|| L oo ((fV „ +a)xR) . 
This bound is similar to the estimate 14.10 I, but has different weights: namely, \w n 



instead of \w 



/+l/2,+ 



- W 



7+1/2,- 



7+1/2,+ 



- W 



7+1/2,- 



. In the first case, (2.9 1 for the subcell states w" +1/2 ± easily yields the 



estimate 

(4.12) w " + i/2,- ~ w 1-i/2,+ = w (u",Vj + i/ 2 ) - w(mJ,i?/-i/2) = 0(l)(Vj+y 2 - Vj-y 2 ) = O(Ax), 



as a consequence of Proposition 4.1 for u^ x and the smoothness property v e W 2/ °°(]R). By contrast, 
the jump in the subcell states \w"^'~_ - w"^'~ + \ cannot be expected to vanish uniformly with Ax, 
since discontinuities may develop within the subcells during the evolution step < |2.16) . One would 
thus expect to control the jump in < |4.9) via a BV estimate, but such an estimate is not available in 
the present framework. This is the reason why we emphasize the sharper inequality | |4.9} over 
the cruder bound 14.10 1. The former will be seen to imply the following estimataM 



1 Such an estimate was first established, for finite difference schemes in several space dimensions, in Coquel and 
LeFloch 1 24 1, who coined the term "weak BV estimate" to denote this bound. The terminology was used in most papers 
on the subject since then. 
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Proposition 4.6 (Entropy dissipation estimate). Let T > be given and Nt be the greatest integer 
smaller than T/At. Then, for any time-independent non-negative test function ip e D(R), the finite 
volume approximation (|2.9|>— (2.30 ' obeys under the CFL condition (2.8 1 the weak BV estimate 



(4.13) 



N T 
n=0 j 



- W 



n+1- I 



fyyAx < 0(1), 



where ipj = l^j-yi + i'j+yz) and ip j+1/2 = ^ f*' +1 \p{x)dx. 

We will use a slightly simpler estimate, depending upon a non-negative test-function if> e 
D(R) such that i/»(x) = 1, |x| < L for some L > 0. Denoting / the largest integer smaller than Lj Ax, 
the estimate 1 4.13) implies 



(4.14) 



n=0 |;'|<; 



4-1 - 

1/2- 



4.3. Convergence arguments. The estimate ( 4.14[ allows us to establish the following (con- 
tinuous in time and space) version of the entropy inequalities. 



Proposition 4.7. Under the CFL condition ( |2.gf , the approximate solutions w& x obey the entropy 
inequality 

J \U(w Ax )d t (p(t,x) + J(w Ax ,v)d x cp + cpL(w Ax ,v)d x v\dtdx + J U(w° Ax )(p(0,x)dx > O(yfKx) 
Jr+xr ^ ' Jr 

for any (smooth) convex entropy pair (U, 3^ in the form ( 2.23 H 2.24 K 

Equipped with the above inequality, we deduce that p~t,x (associated with {itfA*} A Q ) satisfies 

I ((p,,U(.))d t (p(t,x) + (p l ^(.,v))d x cj} + (p,L(.,v)d x v)(p)dtdx + I U(w )(p(0,x)dx > 0. 
Jr+xr^ ' Jr 

In other words, p, is an entropy measure-valued solution of the Cauchy problem 1 1.11) in the 
sense of Definition|3.1| Proving that the initial data \l§ = b Wo with Wq e L°°(R) the initial data of 



the problem l l.ll) is assumed in the strong sense |3.2 1 can be easily deduced from the previous 
analysis by following closely related steps (see for instance |24|). The details are left to the reader. 
By Theorem 3.2 the entropy measure-valued solution p t ,x reduces to a Dirac measure b w (t /X ) 
concentrated on a function w = w(t,x) which coincides with Kruzkov's solution to (1.11 H 1.15) . 
This completes the proof of Theorem 2.1 □ 



Proof of Lemma |43| Under the CFL condition 1 2.8 1, we start from the subcell entropy inequal- 



.,«+!,- 



ities 1 4.4 1 satisfied by the subcell states 
the following entropy inequality centered at Xj + \/ 2 



and W n .\ ' /z . Adding these two inequalities yields 



7+1/2,+ 



M/2 J + U(zv^- + )) 4(U(0^ 1A _) + U(w» +1/2+ )) 

(^■ +1/2/+ ^/+l/2)-^" +1/2/ _ 



Ax\ 



Vj+l/2 



))<o. 



Multiplying this inequality by cp" +l j 2 Ax and summing in space gives 



(4.15) 



Zj Kuizv^-J + Uiw^-^^Ax - Zj l{U(w1 +1/2 J + U(w 
+AtZj(nw1 +V2+ ,v j+U2 ) - nw1 +1/2r _,v j+1/2 ))cp'l +V2 < 0. 



)>f>» „A.v 



;+l/2,+V r i+l/2 



We now deal with the discrete formulation in space and relying on the identities 

H^ +1/2/+ , ^'+1/2) = 5(^+1/2,+' V1/2) + ?(ro? + 3 /2/ _, ^+3/2)) + l(?(w>] +1/2r+ ,Vj +1/2 ) - 3(w'> +m _,v j+3/2 )), 
?(™" +VZ -,v j+1/2 ) = l(^(-w] +1/2 _,Vj +1/2 ) + y(H^_ 1/2/+ ,D;-i/ 2 )) + ^(w'j +1/2 __ r Vj +V2 ) - 3(wn_ 1/2A ,v hl/2 )). 
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We deduce that 



(4.16) 



i 

^ 1 

= ~Yi ^K-lA+'^j-l/z) + 3 r (< + l /2 ,-^/+l/2))(^" + l / 2 " ^-l/ 2 ) 
7 

-£(^"+1/2,-^7+1/2) - n^_ vz+ ,v hV2 )) 2(^+1/2 + ^-i /2 )- 



In view of the definition d2.9| of the subcell states u>" +1/2± an d t ne identity 1 1.14 1 < J{w{u,v),v) = 
Q(u, v), we get on one hand 

l(?( w "-i/ 2 +^7-1/2) + ^"+1/2 _/f;+i/z)) = Q(w ; ",^-i/2) 
(4.17) If 1 ' 

+ 2 J d v Q(Uj,Vj-i /2 + s(v j+y2 - v hV2 ))ds (Vj+i/2 ~ Vj-1/2) = Q(u",Vj-i/ 2 ) + O(Ax). 

Here, we have used the sup-norm estimate in Proposition 1 4.1 1, the smoothness of the mapping 
Q(u, .), together with \Vj + i/ 2 - O/-1/2I = O(Ax), which is a consequence of the property v E W 2 '°°(R) 
and ||23| defining v& x (x). 

On the other hand, by using similar arguments, we obtain 

fiw'j+i/i,-' ^7+1/2) - "-1/2,+' v hVi) = J d v Q(u1,Vj-i/ 2 + s(v j+ i /2 - Vj-i J2 ))ds (v j+ i /2 - Vj-y 2 ) 

= d° v Q(u1,v H / 2 )(v j+1 / 2 - Vf-1/2) + 0((Ax) 2 ). 

(4.18) 

Plugging | |4.17) and | |4.18) in \A.16} then gives us 
Yj {^ wn i+i/2,+' v hV2) ~ ?(™1+ V 2,-'Vj+w))<i>1+V2 



= -^Q( M » / i ?i _ 1/2 ) 



fe/2-^-1/2) 



Ax ~ E i^Vm + (p1-i/ 2 )dvQ(u1,v H/2 ) 



(vj+1/2 - Vj-1/2) 



Ax — i_j 2^7+1/2 ■ rr-wr -r-v Ax 

+0(Ax)||(/)|| wl ,» (R+xK) A>(i",f" +1 ). 
(4.19) 

Finally, routine arguments based on the smoothness of <p and v yield the expected result 



Ax 



(4.20) 



At Y i {^i w 1+i/2,+' v i+V2) - 3(™l + i/ Z -,v j+ i/ 2 ))(p1 +V2 

= — f f ^Q(w A _ Y , v(x))d x cp dtdx + (pd v Q(u& x ,v(x))d x vjdtdx 
+Q(Ax)Af\\cp\\ L ~ (R+xK) x ( t > (t",t" +1 ), 
where, by definition, [ 1.13} -l 2.2 1, one has J(w& x , v(x)) = Q(u& x ,v(x)) and L(w& x ,v(x)) = d v Q(u Ax ,v(x)). 
Proof of Lemma 4.4 We derive (4.8 1 by plugging (2 U(wJ)0J-U(o;^ 1/2 J<£ - +1/2 -U(ri;* 2 + )</>"_ 1/2 ) 



- {^1+1,2,-) + ^1-mM + l M w 1+i/2,-) - u (™U/2,+)){4>1+i/2 - ^/-l/2)' 



with (p n . defined in 1 4. 6 1. In view of 1 2.9 1 (namely, 2w" = iv"_y 2+ + w\y 2 _), Jensen's inequality 
implies 

2 u(wj)cpi - u(w; +1/2 ^ +U2 - u{wi vz+ )r h i, 2 

(4.21) = (2 U(^) - U(^ +1/2/ J - U{w1_ V2+ ))^ + |(lX(^_ 1/2/+ ) - U(^ +1/2 j)(^ +1/2 - ^_ 1/2 ) 

< |(UK- 1/2 , + ) - "K +1/ 2,-))tel/2 - ^"-1/2)' 
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which provides the required bound ( |4.8} in view of the estimate < 4.12} satisfied by the jump 
( w 1+i/2 - ~~ w 'j-i/ 2 +)• Deriving the second estimate ( 4.9 1 is completely analogous, from an identity 
similar to 14.21 \, that is 



(4.22) 



= (2 U(^ +1 ) - U(zv"+lf 2 _) - U(wp y l + ))cp» 



In view of the identity 2w" 
and obtain 



+ 

n+1 _ „,«+! 



w 'j + y 2 - + w j-i'/2+ m 12.17} Lemma 



2.3 



we can apply Lemma 



4.5 



2 U(^ +1 ) - U(^-_) - U(^;- + ) < -joulw^l- - wT-wJ' 



„n+l, 



n+1- 



„n+l- \2 



so that < |4.9} follows from ( 4.22} . This concludes the proof of Lemma 4.4 



Proof of Proposition 4.6 We start from the discrete in time, continuous in space formulation 
4.5} and write 



(4.23) ' r'" +1 r , v 
3 r (w^,y(x))^ v c/)(t,x) + (/)(f,x)£(^ x ,57(x))(9. v z; 

< 0(Ax)Ai||(/)|lwi,»((,,,,, + i)xR)^(f n , f" +1 ), 



in which we plug the decomposition ( |4.7} of the discrete time derivative. Here, we use the discrete 
test function ipj defined in the proposition from any time-independent test function \p e D(R). 
We thus get 

E - UiwfyxjjjAx - J j (?(Kx> v(x))dxW x ) + $(x)£(wlx> vfrWxvJ dtdx 

(4.24) < I E ( 2 U ( w ] +1 Wi ~ U ( w "+i/l-WhV2 ~ Wffii+Wj-w) Ax 
i 

i 

and the estimates 1 4.8 1 and | |4.9} stated in Lemma 4.4 then yield 

E (U{w'] +l ) - U(w?j)^Ax + ^ou E WjSg. - ^lf 2 , + \ 2 4>A x 

i i 
< 0(Ax)A%|| w i,» (IR) + O(Ax) ^ \\Jx-ip\\L«((x HI2 ,x i+V2 ))&x + 0((Ax) 2 ) \\d x xp\\ m ,*, {{x ,_ l/2iX , +1/2)) Ax 

J" ^(a;^., w(x))<9 Y i/>(x) + \l){x)L{w\ x ,v{x))d x v^dtdx. 
The sup-norm estimate in Proposition |4. 1 1 implies the following crude estimate 



(4.25) 



IW!'~(R)/ 



so that 



(4.26) E ( U K +1 ) " Uiw^jhx + ^au E K+fe- " ^ Q{M)\\^\\ w ^ y 
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Summing over all n e [0, Nj] with Nj = [T/ At] (for T > fixed), we get 



f 



U(w Ax (x, T))dx + -o u Y Y \ w "n/2- ~ 

h=0 ;' 



^- + | 2 t/»yAx < f U(w (x))dx + 0(l)T||Vllw^(B)/ 
y Jr 



which gives the desired estimate 14.13 i, choosing for instance the quadratic entropy U(w) = ztr/2 
with ou = 1. This completes the proof of Proposition 4.6 □ 



Proof of Proposition \4.7\ For any e D(R + x R) and in view of its discrete representation 



b", we again consider the continuous in space formulation 1 4.5 i in Lemma 4.3 Plugging in the 



decomposition | |4.7| gives 

Y (U{w"j +l ) - Uiwfy^Ax - J J (?(w n Ax/ v(x))d x( p(t,x) + ^t,x)£(v^Mx))dxvjdtdx 

< 0(Ax)At\\cp\\ w ^ {{t „ jt „ + ^ R) x4t n ,t n+1 ) + 0((Ax) 2 ) Y \\^\\^((t«,t^Mx H/2 ,x i+1/2 )Ax 

i 

+0(Ax) Y 1^2,- - H' 9 ^llL~(( t ", t -)x(, / _ 1/2 ,, /+1/2 ))Ax / 



where we have used the estimates < 4.8 1 and 14.10 1. Summing this inequality over the time indices 
yields 

E (p n+1 -(p n r r / \ 
VuW +1 )-4 -At Ax - I I [3(w n Ax ,v(x))d x ({)(t,x) + (p(t,x)L(w n Ax/ v(x))d x v \dtdx 

n=ax...j At Jr + Jr v ' 

< 0(Ax)\\ct)\\ m ~ (]R+xU) 

+ °(!) Y E - w ]%,M tn > ^ Il5^ll^(( t n, f - I)x (, ; _ 1/2 ^ 1/2 ))A^). 

H>0 j 

Cauchy-Schwarz inequality then allows us to upper bound the last term according to 
Y E { lw hV2,- - ™"%M( tn ' tn+1 )) \\dx<Ph<"MAx 

( E E (Kvi- - ™T-v2,M( tn > t n+i )) 2 AtAxj /2 ( y Y (^A^m^H 

b (t n ,t n+1 )AtAxj 



n>0 j 
< 

n>0 j 

< 0(1) 



H>0 / 



n>0 j 
\l/2 



The weak BV estimate 14.14' implies that L„> Ly \ w 
0(Ax), at most, and routine arguments allow us to conclude. 



n+1- _ - Ji+l- \2 V (in f.n+1 
/+1/2- W ]-V2,+ 



X<p(t n ,f WAx is of order 

□ 



5. Coupling via conservative variable and total variation estimate. This section briefly 
addresses the particular case of the state coupling, namely the case of the coupling condition 
( |1.2| with 6- = 6+ = Id, that is, we impose here w(0-, t) = w(0 + , t) for all t > 0, as considered 
by Godlewski and Raviart [33 1. In this setting, we can derive a BV estimate for the approximate 
solutions w Ax and the convergence follows from Helly's compactness theorem. No smoothness is 
required on the color function v, which can be discontinuous, so that the result in this section holds 
in a stronger norm than the one in Theorem 2.1 but applies to initial data Wq e L°°(]R) nB V(]R), only. 
To avoid technicalities, we restrict attention to the Godunov solver. Note that for the coupling 
zv(0-, t) = w(0+, t), one has u(w, v) = w for all v, and the condition | |5.1} below yields \w\ < ||woIIl°°(R) : 
in view of Proposition |4.1| on e has |MIl°°(r;xR) ^ \\w o||l°°(r), while d w f(w,v) coincides with the 
non-trivial eigenvalue oFJL6}. 
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Theorem 5.1 (Total variation estimate). Assume 8- = 6+ = Id. Let v in L°°(R) and wq in 
L°°(]R) n BV(]R) be the initial data. For any small e > (e < 1/2) and under the strengthened CFL 
condition 

1 

(5.1) max \d w f(w,v)\<--e, 

|i»|<||i» || L oo (R) Z 

H<1 



the finite volume method {2.9 h |2~T0 > is total variation diminishing and, in particular, TV(w Ax (t", ■)) < 
TV(w Q ). 

Proof. The case of the state coupling comes with y_ = y+ = Id, therefore the reconstruction 
step boils down to w"_ 1 / 2+ = w> j + i/2- = w 'j' nam ely the discrete solution is kept unchanged and 
no jump is created at this step. Consider the Riemann problem at the interface X/+i/2 (with Cy+1/2 
constant) and denote w(-,w n .,w n . + ^) its self-similar solution, which is known to be monotone and 

satisfies TV(w(.,w" ,w" +1 )) = \w" +1 - w% Define w Ax {t n+l ~ ,x) by glueing together non-interacting 



neighboring Riemann solutions. Under the CFL condition 1 5.1 1, we find 

TV(w Ax (t n+l -, •)) = Lj TV lXj+eAx /2 Aj+1 -eA x /2](lV Ax (t" +1 -, •)) + Lj TV (xreAx /2, Xi+ eA x /2)(lVA x (t n+1 -, •)) 
= LjTViwi^u'j,™^)) + IZjTVty-eAx/lXj+eAx/ZiiWAxit",-)), 

but the second term vanishes since w"_ lj2+ = w" +1/2 _. Therefore, TV(zv Ax (t" +1 ~ , ■)) = £y \w" +l - w"\ 
= TV(w Ax (t n ,■)). Denote 7 Ax the operator defined (for x e (x y _i/2,x ;+ i/2)) by w Ax (t n+l ,x) = 
T 'ax(u> 'Ax(t n+1 ~ ' , x)) = ^ Ix 'tfz W ^' l+1 ~' y) dy- This averaging operator Tax is total variation dimin- 
ishing, and so TV{w Ax {t n+1 , ■)) < TV{w Ax {t n , ■)). a 
6. Numerical experiments. 

6.1. Formulation of the problem. In this section, coupling problems are numerically studied. 
Resonance in the infinitely sharp regime is closely investigated for smooth and discontinuous 
coupled solutions. The color function modeling the thick interface is chosen to be v n {x) = 
(erf(x/?]) + l)/2 (x e IR) for some fixed rj. Here, x i-» erf(x) denotes the classical error function. 
Obviously v n takes value in the expected range (0, 1) and belongs to W 2 '°°(]R) (in fact to e°°(]R)). 



Observe that x i-» c,,(x) - 1 is an even function so that the coupling formula 1 1.7 1 is, in some sense, 
"symmetric" in the left- and right-hand partial differential equation. In order to investigate the 
dependency of the discrete solutions upon the regularization, we consider the (two parameters) 
function v, v q(x) = (erf(x/n + Q + l)/2 (x e IR) with 77 > and CeK. Here, 1] clearly monitors the 
thickness of the handshake coupling zone while C acts roughly speaking as a shifting perturbation 
breaking the symmetry of the treatment of the left and right problems. All the calculations are 
performed over the bounded domain [—1,1] and Neumann boundary conditions are used. 

6.2. Non-uniqueness for a resonant coupling problem. We consider the closure relations 

f~(w) - y, f + (w) - ^ W+ }' , 6-{w) = 6+{w) = w, and the initial data is Wo(x) — Wi — —\ for x < 0, 
and w r = 3/2 for x > 0. The corresponding Riemann solution in the regime of an infinitely thin 



interface exhibits a resonance phenomenon which is depicted in Figure 6.1 Observe that in the 
present problem, the sonic point associated to the right flux / + (respectively to the left flux f~) is 
u s + = -1 (resp. U s _ = 0). The ordering of the left- and right-hand sonic states: Mi > w+implies a 
resonance phenomenon in outgoing waves from the coupling interface. This is a sign of failure of 
uniqueness in solutions of the proposed coupled Riemann problem. We refer the reader to [17 1, 
the references therein, and to Figure [6T| 

In the regime of a regularized thick coupling interface, uniqueness of the solution of the 
Cauchy coupled problem is restored. A given regularization actually selects a given intermediate 
state so as to define a unique solution. Our aim here is to provide numerical evidences supporting 
this claim. We use the two-parameters family of regularizations c, ; -q described above, for various 
thickness i] and shift Q. Discrete solutions have been computed under the CFL condition ( |2.8) 
using respectively 100 and 1000 grid points. The regularization parameters are respectively set to 
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the constant value t] = 5.10 -3 and to three different values concerning Q. respectively -0.5,0 and 
0.5. The extreme values select the left- and right-hand problems depending of sgn(C) (negative 
values select /~), while the intermediate value keeps a symmetry within the two problems. The 
numerical results are displayed in Figures 6.2 and 6.3 Two facts must be highlighted: the mesh 
refinement must be fine enough so as to capture a constant intermediate state, then and clearly 
the resulting value of the intermediate state it>* is highly sensitive with respect to the choice of the 
shifting parameter C- In Figure 6.4 we investigate the dependence with respect to the thickness 
parameter t], while the shifting parameter C is kept fixed at the constant value 0.5. Two values of 
1] are considered: 1] = 0.01 and rj = 0.001. The resulting w profiles obtained using 5000 grid points 
turn significantly less sensitive to the choice of the thickness parameter rj than to the choice of the 
shifting parameter £. 




W+ e [w e , w r ] 

Figure 6.1: A one-parameter family of smooth solutions in a resonant situation 



interface - zeta= 0.0 - 
interface - zeta= 0.5 
interface - zeta=-0.5 ■ 



solution - zeta - 0.0 - 
solution - zeta = 0.5 
solution - zeta =-0.5 ■ 
u at t=0.0 




(a) Numerical color functions (b) Numerical solutions 

Figure 6.2: Failure of uniqueness in a resonant situation - N — 100 



6.3. Another example of non-uniqueness. This last test case is based on the closure relations 
in Section 6.2 but the initial data Wq{x) = W[ = lforx < Oandzf, = -2forx > 0. Such a choice again 
results in a resonance phenomenon involving multiple solutions. Three different discontinuous 
solutions are available in the regime of coupling interface with zero thickness. We refer the reader 
to IIT61 IT7I for more details. These discontinuous solutions, depicted in Figure 6.5 respectively 
coincide with a left moving shock, namely a discontinuity satisfying the Rankine-Hugoniot jump 
relation for f~, a standing discontinuity, i.e. a discontinuity with zero speed, and a right moving 
shock, that is a discontinuity satisfying the Rankine-Hugoniot jump relation for / + . We show 
hereafter that these three different solutions are actually stable, in the sense that we can capture 
each of them numerically by choosing suitable regularization parameters in the thick coupling 
interface v^c- We choose as previously a fixed thickness rj = 5.10 -3 and we make the shifting 
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interface - zeta = D.O - 
interface - zeta = 0.5 
interface - zeta =-D. 5 ■ 



solution - zeta = 0.0 

solution - zeta - 0.5 

solution - zeta =-0.5 

u at t=0.0 








y 









0.1 -1 



(a) Numerical color functions (b) Numerical solutions 

Figure 6.3: Failure of uniqueness in a resonant situation - N = 1000 




-0.06 -0.04 -0.02 0.02 0.04 0.06 -1 -0.5 0.5 1 



(a) Numerical color functions (b) Numerical solutions 

Figure 6.4: Sensitivity with respect to the thickness parameter 



parameter C to vary from the negative value -0.5 (selecting the left-hand equation) to the positive 
one 0.5 (selecting the right-hand equation) and the vanishing value C = 0, keeping the symmetry 
in between the two. As heuristically expected, the results displayed in Figure 6.6 show that the 
negative value of C captures the left moving shock, the right-hand state selects the right moving 
shock, while restores the standing discontinuity. In particular, multiple solutions may be stable. 
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(a) A left shock (b) A standing discontinuity (c) A right shock 

Figure 6.5: Several discontinuous solutions in a resonant situation 
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solution - zeta = 0.0 

solution - zeta = 0.5 

solution - zeta =-0.5 



-0.1 -0.05 0.05 0.1 -1 -0.5 0.5 1 

(a) Numerical color functions (b) Numerical solutions 

Figure 6.6: Failure of uniqueness in a resonant situation with discontinuous solutions - N = 1000 
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